Coded spectral imager

ABSTRACT

Embodiments herein provide for imaging objects. In one embodiment, a spectral imaging system includes an optical element configured to receive electromagnetic energy of a two-dimensional scene and a filter configured to provide a plurality of spectral filter profiles. The filter also transmits multiple spectral wavebands of the electromagnetic energy substantially simultaneously through at least one of the spectral profiles. The spectral imaging system also includes a detector configured to measure intensities of the multiple spectral wavebands, and a processor configured to generate a spectral image of the scene based on the measured intensities.

CROSS REFERENCE TO RELATED APPLICATIONS

This patent application is a non-provisional patent application claimingpriority to, and thus the benefit of an earlier filing date from, U.S.Provisional Patent Application No. 62/042,570 (filed Aug. 27, 2014), theentire contents of which are hereby incorporated by reference.

FIELD OF THE INVENTION

The invention relates to spectral analysis of images.

BACKGROUND

A scene viewed by a hyperspectral imaging system is often displayed as athree-dimensional datacube, where two dimensions represent the spatialdomain (x, y) and one dimension represents the spectral domain (λ).Spatial observations (x, y) allow a person to observe an image when highcontrast is available. However, during conditions of low contrast, suchas fog, smoke, camouflage, and/or darkness, spectral signatures helpidentify otherwise unobservable objects. Hyperspectral imagers arecapable of collecting and processing objects within a scene by dividingthe whole spectrum into tens or even hundreds of bands and thus canobtain high resolution datacubes useful in a wide range of applicationssuch as mining, agriculture, chemical detection, and militarysurveillance. Conventional hyperspectral devices face a tradeoff betweenspectral resolution and the signal to noise ratio (SNR) of the estimatedspectrum. Thus, there is an ongoing need for imaging systems thatmaximize the signal to noise ratio (SNR) of the estimated spectrumwithout sacrificing capability of generating datacubes with highspectral resolution.

SUMMARY

Systems and methods presented herein provide for imaging an object. Inone embodiment, a spectral imaging system includes an optical elementconfigured to receive electromagnetic energy of a two-dimensional sceneand a filter configured to provide a plurality of spectral filterprofiles. The filter also transmits multiple spectral wavebands of theelectromagnetic energy substantially simultaneously (i.e., at or aboutthe same time) through at least one of the spectral profiles. Thespectral imaging system also includes a detector configured to measureintensities of the multiple spectral wavebands, and a processorconfigured to generate a spectral image of the scene based on themeasured intensities.

The various embodiments disclosed herein may be implemented in a varietyof ways as a matter of design choice. For example, some embodimentsherein are implemented in hardware whereas other embodiments may includeprocesses that are operable to implement and/or operate the hardware.Other exemplary embodiments, including software and firmware, aredescribed below.

BRIEF DESCRIPTION OF THE FIGURES

Some embodiments of the present invention are now described, by way ofexample only, and with reference to the accompanying drawings. The samereference number represents the same element or the same type of elementon all drawings.

FIG. 1 is a block diagram of an imaging system in an exemplaryembodiment.

FIG. 2 is a flowchart illustrating an exemplary process of the imagingsystem of FIG. 1.

FIG. 3 is a flowchart illustrating an exemplary process of the imagingsystem of FIG. 1.

FIG. 4 is a block diagram of an imaging system in another exemplaryembodiment.

FIG. 5 illustrates transmission characteristics of multiple-wavebandspectral filters in an exemplary embodiment.

FIG. 6 is a block diagram of an imaging system in another exemplaryembodiment.

FIG. 7 is a block diagram of an imaging system in another exemplaryembodiment.

FIG. 8 illustrates a QRS sequence and the result of convolving the QRSsequence with its deconvolution kernel.

FIG. 9 illustrates a graph of the gain ratio for a variety of values.

FIG. 10 is a graph that illustrates the SNR gain when exploiting theanti-correlated nature of the noise.

FIG. 11 illustrates a graph of the SNR gain as a function of χ.

FIG. 12 illustrates an exemplary sequence of images taken with the samefilter which is circularly shifted for each image.

FIG. 13 illustrates an exemplary panchromatic sequence of images takenof a scene which contains four circular sources each with a uniquespectral profile.

FIG. 14 illustrates actual spectrum for each source of FIG. 13 in anexemplary embodiment.

FIG. 15 illustrates a high resolution reconstruction of the spectrumfrom the panchromatic images of FIG. 13 in an exemplary embodiment.

FIG. 16 illustrates a reconstruction of the spectrum from thepanchromatic images of FIG. 13 in an exemplary embodiment.

FIG. 17 illustrates a noisy image in an exemplary embodiment.

FIG. 18 illustrates a comparison of spectrums for a pixel in the noisyimage of FIG. 17 in an exemplary embodiment.

FIG. 19 illustrates a comparison of spectrums for a pixel in the noisyimage of FIG. 17 in another exemplary embodiment.

FIG. 20 illustrates images obtained from simulated measurements of ahyperspectral datacube in an exemplary embodiment.

FIG. 21 illustrates a comparison of the spectrum for the central pixelof FIG. 20.

FIG. 22 is a block diagram of an exemplary computing system in which acomputer readable medium provides instructions for performing methodsherein.

DETAILED DESCRIPTION OF THE FIGURES

The figures and the following description illustrate specific exemplaryembodiments of the invention. It will thus be appreciated that thoseskilled in the art will be able to devise various arrangements that,although not explicitly described or shown herein, embody the principlesof the invention and are included within the scope of the invention.Furthermore, any examples described herein are intended to aid inunderstanding the principles of the invention and are to be construed asbeing without limitation to such specifically recited examples andconditions. As a result, the invention is not limited to the specificembodiments or examples described below.

FIG. 1 is a block diagram of an exemplary imaging system 100. Imagingsystem 100 is capable of generating hyperspectral images of a scene 102,where each hyperspectral image includes a full spectral profile (e.g.,light intensity as a function of wavelength) of an individual pixel of atwo-dimensional image. A series of hyperspectral images may be organizedinto a three-dimensional datacube, where two dimensions represent thespatial domain (x, y) and one dimension represents the spectral domain(λ) of scene 102.

Imaging system 100 includes an optical element 110, a filter 120, adetector 140, and a processor 150. Optical element 110 is anycombination of devices operable to propagate electromagnetic energy(e.g., photons) from an object or scene 102 through imaging system 100.Such devices may include focusing lenses, relay lenses, collimatinglenses, dispersers, gratings, reflectors, etc.

Filter 120 is any system, component, or device operable to transmit intwo or more wavebands simultaneously, or substantially simultaneously,and to attenuate in at least one other waveband. In one embodiment,filter 120 comprises multiple spectral profiles 122-126, any of whichare operable to pass multiple distinct wavebands, or wavelengthspectrums. For simplicity, FIG. 1 shows filter 120 implemented withthree spectral profiles 122-126, each transmitting differentcombinations of two out of three spectral bands, marked λ1, λ2, and λ3.However, alternative configurations, number of spectral profiles122-126, and/or number/combination of transmitted wavebands through eachspectral profile 122-126 are possible in accordance with the particularspectral range of interest. Additionally, as will be discussed infurther detail below, imaging system 100 may implement one or multiplefilters in alternative configurations.

Detector 140 is any system, component, or device operable to receive themultiple wavebands of electromagnetic energy transmitted through filter120 and to generate one or more electrical signals that indicateelectromagnetic intensity as a function of wavelength. As used herein,detector 140 may refer to a single detector element, or pixel, within atwo-dimensional detector array of detector elements. Alternatively,detector 140 may refer to the detection of a larger area or entiretwo-dimensional image of scene 102. Detector 140 may comprise, forexample, a charge-coupled device (CCD) sensor, a silicon basedcomplementary metal-oxide-semiconductor (CMOS) sensor, or other types ofsensors and sensor elements.

Processor 150 is any system, component, or device operable toreconstruct a hyperspectral image of scene 102 using the intensityvalues of detector 140. Thus, processor 150 may be communicativelycoupled with detector 140 to receive the intensity values output bydetector 140. Processor 150 may use digital processing techniques todemodulate the measured intensities into separate, individual wavebandsignals. Thus, though multiple lines of radiation are focused onto asingle detector or detector element, processor 150 may separatelyanalyze the characteristics of each line of radiation for use inanalyzing the composition of scene 102.

Previous imaging systems obtain hyperspectral images at high spectralresolution by using a narrow passband, or narrow spectrometer inputslit, and scanning one wavelength band sequentially over time. Sincespectral information in such imagers is not acquired at the same time,these imagers are not able to benefit from the improvement in signal tonoise ratio (SNR) that multiplexing provides. The SNR is further reducedin these systems because the passband reduces the throughput ofelectromagnetic energy through the imager. For example, if an incidentsignal with a uniform spectrum is split into 100 equally separatedwavebands, the photon flux into each waveband will be 1% of the photonflux from the entire signal. Consequently, the SNR for a detectormeasuring a single waveband would be 10 times less than the signal tonoise ratio for a detector measuring the entire photon flux associatedwith all the wavebands.

Other spectral imaging systems use a multiplexed focal plane array wherean image is passed through an interferometric assembly. This technique,commonly referred to as Fourier transform spectrometry, enablesmultiplexing of a two-dimensional spatial image without limiting thelight throughput with a narrow slit. However, these systems record thespectral dimension with precise mirror positioning which makes suchsystems susceptible to vibrations and therefore unsuitable for imagingsystems in airplanes, unmanned aerial vehicles, and other applications.

Imaging system 100 advantageously enables a combination of highthroughput of radiation and simultaneous multiband or multiplexwavelength scanning. The higher throughput arises from the fact thatfilter 120 of imaging system 100 allows a larger transmissive area thanthe slits or passband filters, thus enabling higher throughput ofradiation, increased signal measurement, and larger SNR ratios. Further,filter 120 and detector 140 allow imaging system 100 to collectelectromagnetic radiation from a plurality of spectral bands insimultaneous fashion, therefor providing imaging system 100 with amultiplex advantage without complication of components susceptible tovibrations.

Imaging system 100 also allows for collection of two-dimensional images,wide-field spectral images over time. Additionally, as discussed ingreater detail below, design characteristics of imaging system 100 allowprocessing functions to improve post-processing of the resultingdatacube. Imaging system 100 is thus useful for a variety of imagingapplications such as persistent surveillance where high spectralresolution is desired (e.g., hyperspectral applications) and/or wherethere is a limit to the amount of electromagnetic energy available(e.g., low light conditions).

FIG. 2 illustrates a flowchart of an exemplary process 200 of imagingsystem 100. In step, 202 optical element 110 receives electromagneticenergy of a two-dimensional scene (e.g., scene 102). Optical element 110may focus a two-dimensional image of scene 102 in the plane of filter120. In step 204, filter 120 transmits multiple spectral wavebands ofthe electromagnetic energy substantially simultaneously through at leastone of a plurality of spectral profiles 122-126. In step 206, detector140 measures intensities of the multiple spectral wavebands transmittedthrough the spectral profile. In step 208, processor 150 generates aspectral image of the scene based on the measured intensities.

FIG. 3 illustrates a flowchart of another exemplary process 300 ofimaging system 100. FIG. 3 describes processes which may be performed byprocessor 150 during or after generation of a spectral image (e.g., step208 of process 200). In step 302, processor 150 initiates analysis of adigital output signal from detector 140. In step 304, processor 150selects a sensor element of detector 140, or pixel, for processing.

In step 306, processor 150 generates the spectrum of a pixel in thespectral image by demodulating the encoded spectral information on thatsensor element according to the properties of imaging system 100.Processor 150 is operable to reconstruct underlying information of scene102 using algorithms which decode the multiplexed information receivedat detector 140 using properties of filter 120 and/or other elements inimaging system 100. In that regard, processor 150 may receive, store, orbe programmed with information that models imaging system 100 to enableprocessor 150 to generate spatial and/or spectral information of a pixelbased on the intensity and position of the selected sensor element.

In step 308, processor 150 may optionally perform post-processing on thehyperspectral image to reduce errors in the image. Additional detail ofsuch post-processing will be described in greater detail below. Then, instep 310, processor 150 generates a datacube of scene 102 whichcomprises an array of spectral images (e.g., obtained by iterativelyperforming process 200 and/or steps 302-308).

It will be appreciated that one or more steps of processes 200/300 maybe performed on imaging system 100 having alternative configurations.FIG. 4 is a block diagram of imaging system 100 with multiple staticfilters 120-1 and 120-2 and multiple cameras 440-444 or detectors in anexemplary embodiment. Static filters 120-1 and 120-2 each transmit aportion of the spectrum (e.g., light 450) and reflect another portion ofthe spectrum. That is, the exemplary imaging system 100 of FIG. 4implements a series of static multiple-waveband filters thattransmit/reflect a different portion of the spectrum from one another.

FIG. 5 illustrates transmission characteristics of differentmultiple-waveband spectral filters in an exemplary embodiment. Plot 510shows the transmission characteristics of a filter that transmitsfrequencies of green (G) and blue (B) light and attenuates frequenciesof red (R) light. Plot 520 shows the transmission characteristics of afilter that transmits RG frequencies and attenuates B frequencies. Plot530 shows the transmission characteristics of a filter that transmits RGfrequencies and attenuates G frequencies.

Returning to the example imaging system 100 of FIG. 4, filter 120-1transmits spectrum T1 to camera 440 and reflects spectrum R1 towardfilter 120-2. That is, filter 120-1 transmits 50% of both the R and Bwaveband, while reflecting all of the G waveband and 50% of both the Rand B waveband. Camera 440 receives light combined in the R and Bwavebands. Filter 120-2 transmits spectrum R1×T2 to camera 444 andreflects spectrum R1×R2 to camera 442. That is, filter 120-2 receiveslight that has been reflected from filter 120-1 and transmits all of theincident R waveband and 50% of the incident G waveband. Light in the Gand B wavebands are reflected from filter 120-2 to result in an image atcamera 442. The combination of the reflecting filter properties offilter 120-1 and the transmitting properties of filter 120-2 result inan image at camera 444 comprised of the R and G wavebands.

In this example, imaging system 100 does not substantially discard lightwithin any of the sensed wavebands, but routs light with differentwaveband combinations to different image planes and/or in different ororthogonal directions. Cameras 440-444 may be positioned to coincide atsubstantially identical virtual image planes associated with opticalelement 110. Pixel coordinates from each camera 440-444 may beregistered relative to each other so that a sensed signal magnitude fromeach camera may be associated with a single location at the imaged scene102. Thus, cameras 440-444 and filters 120 may be arranged to align suchthat different combinations of the spectrum from the same location in animaged scene are incident on different sensor elements (e.g., pixels) orcamera 440-444 at the same time or substantially the same time.Processor 150 (not shown in FIG. 4) may implement a processing model ofimaging system 100 to determine, for example, the red, green, and bluespectral components of the image.

Alternatively, for distant scenes, imaging system 100 may directmultiple separate cameras at the same scene. FIG. 6 is a block diagramof imaging system 100 in another exemplary embodiment. Here, imagingsystem 100 includes a RB camera 640 that detects images in the red andblue spectrum, a BG camera 642 that detects images in the blue and greenspectrum, and a GR camera 644 that detects images in the green and redspectrum. In the example imaging system 100 of FIG. 6, each camera640-644 is fitted with its own multi-waveband spectral profile filterand the cameras 640-644 are registered so that light from the samelocation in the scene that is detected by each of the cameras may becompared by processor 150.

FIG. 7 is a block diagram of imaging system 100 in yet another exemplaryembodiment. As shown in FIG. 7, imaging system 100 may additionallyinclude coded aperture 130 and/or optical element 112. Imaging system100 may also include one or more dynamic filters 120 controllable by adriver such as controller 160.

Coded aperture 130 is any system, component, or device operable toprovide aperture coding to multiplex images. Coded aperture 130 maycomprise a coded matrix in the form of a two-dimensional pattern ofsections, where each section may block or reflect some or all of theimpinging electromagnetic energy to accomplish a transform of thespectral data within each of the transmitted wavebands. Coded aperture130, either alone or in combination with other devices such as opticalelement 112, may generate a dispersed mask image comprising multipleimages of the coded aperture corresponding to the different transmittedwavebands. Thus, coded aperture 130 and/or optical element 112 mayinclude a dispersive element (not shown) to induce awavelength-dependent spatial shift of the multiple wavebands so that thetransmitted wavebands may be spatially dispersed onto detector 140.

Detector 140 may comprise a two-dimensional detector array of detectorelements, where each column in the array includes linear combinations ofcolumns of coded aperture 130 that depend on the spectrum of the source.A single detector element, or pixel, may measure multiple spatiallyoverlapped or simultaneously transmitted wavebands. Processor 150 mayretrieve the spectrum of each pixel of a hyperspectral image usingmathematical properties of filter 120, coded aperture 130, and/ordetector 140. Thus, the algorithms or processing of decoding amultiplexed waveband at a single detector element may be dependent onthe particular design of imaging system 100. Processor 150 may perform atransformation of the data matrix to produce a reconstructed spectralimage that mathematically represents the optical spectrum of source 102.

It will be appreciated that imaging system 100 may incorporate anycombination of the exemplary embodiments described and shown herein. Forinstance, imaging system 100 may incorporate a static filter approach, adynamic filter approach, or some combination thereof. Additionally,alternative arrangements, combinations, and number of filters 120,profiles 122-126, and/or cameras 440-444 are possible as well asalternative ranges of electromagnetic energy. For example, the conceptsdescribed herein may apply to detection of other filterable energyfields such as acoustic energy, and particle detection, as well as otherelectro-magnetic energy fields such as radio waves, terahertz radiation,x-rays, and millimeter waves. For particle detection, filteredparameters may include mass, spin, electrical or magnetic moments, orenergy. Furthermore, though imaging system 100 of the exemplaryembodiments may be described and shown herein as having a particularphysical order and location of components, it will be appreciated thatdesigns with alternative orders, respective locations of elements,and/or with a fewer or greater number of components are possible.

EXAMPLE

For the purposes of discussion and analysis that follows, consider thatfilter 120 of FIG. 1 comprises an acousto optical tunable filter (AOTF).In an AOTF, a radio frequency (RF) input applied to a transducercontrols the transmitted radiation intensity level of filter 120.Imaging system 100 may thus include or be communicatively coupled with acontroller (e.g., controller 160) that varies the RF frequencycorresponding to the desired wavelength range to acquire a completespectrum for analysis.

Filter 120 may transmit multiple, discrete wavebands by driving the AOTFat multiple spaced frequencies. Alternatively or additionally, eachspectral profile 122-126 of filter 120 may be segmented into discretezones which are independently controlled by one or more AOTF transducers(e.g., controller(s) 160). Thus, imaging system 100 may select multiplewavebands (or colors) simultaneously.

The image formation process far the image i^(th) image in a sequence ofimages can be described quite generally by the linear superpositionintegral

I _(i)(u,v)=∫_(t) _(i) ^(t) ^(i+1) ∫_(λ) ₀ ^(λ) ¹ ∫_(η) ₀ ^(η) ¹ ∫_(ξ) ₀^(ξ) ¹ PSF(u−ξ,v−η,λ)U ₀(ξ,η,λ,t)×S(λ,t)T(t)dξdηdλdt+n_(i)(u,v),  Equation (1):

where I_(i) (u, v) is the image formed by the imaging system at focalplane coordinates (u, v) at time interval t_(i)≦t≦t_(i+1),PSF(u−ξ,v−η,λ) is the system point spread function (PSF), U₀(ξ, η, λ, t)is the object being imaged at the object coordinates (ξ, η) at a givenwavelength λ and which may depend on time for non-static scenes, S(λ, t)is the spectral response of the detector which may change in time (e.g.,filter 120 comprises an AOTF and the spectral response is tunable), andT(t) is the temporal response of the detector. Here, we may take T(t)=1while the shutter is open and T(t)=0 when the shutter is closed, thus wecan remove the T(t) term from the integral and simply use the limits ofthe time integration to account for an open and closed shutter. Thoughthe wavebands may be discussed herein with respect to models of evenlyspaced spectral bands for convenience of discussion, it will beappreciated that imaging system 100 may select wavebands ofalternatively spaced spectral bands.

The limits of the integration extend over the field of view of thecamera for ε and η, over the spectral bandpass of the detector (e.g.,450-800 nm for a visible band camera), and the time integral extendsover the integration time of the camera T=t_(i+1)−t_(i). The variablen_(i)(u, v) denotes the noise in the i^(th) image such as shot noise,dark current noise, read out noise, etc. Here, assume that the noise isadditive and that the PSF in Eq. (1) is written in the space-invariantform for simplicity. To further simplify, the scene U_(o) may consideredindependent of time. Additionally, the spectral filter may be assumedconstant during the time period that the shutter is open (i.e., whent_(i)≦t≦t_(i+1)). The spectral features of Eq. 1 may be defined as

U′ ₀(u,v,λ)=∫_(t) _(i) ^(t) ^(i+1) ∫_(λ) ₀ ^(λ) ¹ ∫_(η) ₀ ^(η) ¹ ∫_(ξ) ₀^(ξ) ¹ PSF(u−ξ,v−η,λ)U ₀(ξ,η,λ)dξdηdλdt  Equation (2):

This definition along with the simplifications just mentioned allows Eq.(1) to be rewritten as

I _(i)(u,v)=∫_(t) _(i) ^(t) ^(i+1) ∫_(η) ₀ ^(η) ¹ U′ ₀(u,v,λ)S(λ)dλ+n_(i)(u,v),  Equation (3):

where S_(i)(λ)=S(λ, t_(i)) is the spectral filter used for the i^(th)image.

Previous systems implement a narrow band spectral response such thatS_(i)(λ)≈δ(λ−λ_(i)). In that case, the imager may take N snap shots fordiffering values of λ_(i) to build a discrete approximation of U′₀(u, v,λ) with N spectral bands. A standard color camera operates in such a waywith N=3 and λ_(i) is chosen to be red, green, and blue (RGB). Becausethe RGB spectral bands are well separated one can have fairly broadbandfilters which allow a sufficient amount of light to the detector withoutmuch overlap between the spectral transmission lines. However if highspectral resolution is desired (e.g. such as for hyperspectral imagingapplications), the spectral filters may be narrow so as to not overlapand not cause a lot of spectral blurring. This restricts the amount oflight reaching the sensor and lowers the signal to noise (SNR).

Here, imaging system 100 captures multiples images I_(i) (u, v) with atunable filter 120 in order to determine the datacube (two space and onespectral dimension) U′_(o)(u, v, λ) defined in Eq. (2). A tunable filtermay be modeled as

S _(i)(λ)=Σ_(j=0) ^(N-1) q _(ij)δ(λ−λ_(j))*g(λ),  Equation (4):

where q_(ij) is an M×N (e.g., M=N) matrix consisting only of 1s and 0simplying whether passband is present or not present respectively,δ(λ−λ_(j))*g(λ) is a delta function convolved with a Gaussian or similarshaped function, g(λ), that models the passband of the spectral filter,and λ_(j) is the center frequencies of each passband in the filter.Substituting the expression into Eq. (3) gives us

I _(i)(u,v)=Σ_(j=0) ^(N-1) q _(ij) U″ _(o)(u,v,λ _(j))+n_(i)(u,v),  Equation (5):

where U″_(o)(u, v, λ)=U′_(o)(u, v, λ)*g(λ), is the original datacubespectral convolved with the filter passband function g(λ). Thus, Eq. (5)is a matrix equation for the spectrum at a given spatial coordinate (u,v) on a two-dimensional array of detector 140. Here, detector 140 maycomprise a charge-coupled device (CCD) sensor, a silicon basedcomplementary metal-oxide-semiconductor (CMOS) sensor, or other types ofsensors. Filter 120 may implement a pattern that defines each column ofan invertible coding matrix. To solve for the spectral datacube, imagingsystem 100 takes a sequence of images with different filters such thatq_(ij) is invertible. Because there are multiple passbands transmittedthrough filter 120 (and coded aperture 130 if implemented), the samespectral line is measured multiple times at an element of detector 140.Processor 150 may invert the matrix is inverted to solve for thespectral datacube. In the process of inverting the matrix, each of themultiple measurements add together for a √{square root over (N)}improvement in SNR.

Suppose, in one embodiment, that a filter is chosen for the imageI_(i+j)(u, v) to be a circularly shifted version of the filter I_(i).Then the matrix q_(ij) is a circulant matrix such thatq_(ij)=q_(j-i (mod N)). The main advantage of using circulant matrices,is that their inverse is known, and may be quickly calculated usingFourier transforms. Thus, processor 150 may use discrete Fouriertransform processing to demodulate the signal to determine, in realtime, the intensity of each line of radiation reflected onto an elementof detector 140.

The solution estimate for U″o(u, v, λ_(j)) in Eq. (5) is

$\begin{matrix}{{U_{est}\left( {u,v,\lambda_{I}} \right)} = {\frac{1}{N}{\sum\limits_{j = 0}^{N - 1}{\sum\limits_{k = 0}^{N - 1}{\frac{I_{k}\left( {u,v} \right)}{\psi_{j}}e^{{- 2}\pi \; {i{({k - l})}}{j/N}}}}}}} & {{Equation}\mspace{14mu} (6)}\end{matrix}$

where Equation (7): ψ_(j)Σ_(m=0) ^(N-1)q_(m)e^(2πimj) ^(N)

is the inverse discrete Fourier transform of one row of the matrix,which are also the eigenvalues of the circulant matrix. The estimator isobtained by inverting the matrix q in Eq. (5). The estimator would beexact in the absence of noise.

Suppose, in another embodiment, that the filter may be further definedsuch that not only is the matrix q_(ij) of coded aperture 130 circulantbut it is also a quadratic residue sequence (QRS). For quadraticresidues there exists a basic function

$r_{i} = {\frac{2}{N + 1}\left( {{2q_{i}} - 1} \right)}$

such that

Σ_(i=0) ^(N-1) q _(i) r _(i-j)=δ_(ij).  Equation (8):

Processor 150 may use this identity to estimate the datacube U″_(o)(u,v, λ_(j)) in Eq. (5) by the operation

U _(est)(u,v,λ _(j))=Σ_(i=0) ^(N-1) I _(i)(u,v)r _(i-j) =U″ ₀(u,v,λ_(j))+Σ_(i=0) ^(N-1) n _(i)(u,v)r _(i-j),  Equation (9):

Quadratic residue sequences are discrete sequences that are binary,positive functions containing only ones and zeros of roughly equalnumber. Thus, a QRS may describe the transmission profile ofapproximately 50% transmissivity. However, because QRS have the propertyshown in Eq. (9), they may be easily inverted so that they act as asingle filter with the spectral resolution of the smallest open element.Additionally, the matrix is well conditioned and the noise is notseverely amplified.

FIG. 8 illustrates a QRS sequence and the result of convolving the QRSsequence with its deconvolution kernel. Plot 810 shows a graph of alength 31 QRS filter convolved with a 5 nm width Gaussian. The QRSsequence is expanded to a discrete grid 310 units long, and convolvedwith a Gaussian with a standard deviation of 5 nm. This mimics thefilter properties of an Acousto Optical Tunable Filter. Plot 820 showsthe result of convolving the length 31 QRS filter function with itsdeconvolution kernel r_(i), as described in Eq. (8). The result is theGaussian filter function of standard deviation of 5 nm, which is thesame function that was convolved with the original sequence.

To solve for the spectral datacube when using tunable filter 120 andcirculant QRS coded aperture 130, processor inverts the matrix qij. Ifthe matrix is ill-conditioned then the act of inverting Eq. (5) mayamplify the noise and eliminate the SNR gain we expect from takingmultiple images of the same spectral band. To see the effects that noisehas on the signal estimation process we may substitute Eq. (5) into Eq.(6). The estimate for the spectral datacube in the presence of noise is

$\begin{matrix}{{U_{est}\left( {u,v,{l\; \lambda_{m}}} \right)} = {{U_{o}^{''}\left( {u,v,\lambda_{j}} \right)} + {\frac{1}{N}{\sum\limits_{j = 0}^{N - 1}{\sum\limits_{k = 0}^{N - 1}{\psi_{j}^{- 1}e^{{- 2}\pi \; {i{({k - l})}}{j/N}}{{n_{k}\left( {u,v} \right)}.}}}}}}} & {{Equation}\mspace{14mu} (10)}\end{matrix}$

Assuming that all the noise terms are independent Gaussian processeswith mean zero allows us to calculate the total variance of the signalestimate in Eq. (10) as:

$\begin{matrix}{{\sigma_{est}^{2} = {\frac{1}{N^{2}}{\sum\limits_{j = 0}^{N - 1}{\sum\limits_{k = 0}^{N - 1}{{\psi_{j}}^{- 2}\sigma_{k}^{2}}}}}},} & {{Equation}\mspace{14mu} (11)}\end{matrix}$

where σ² _(k)=E[n_(k)(u,v)²] and E[ . . . ] implies the expectationvalue obtained after a series of many measurements. The functionaldependence on (u,v) is dropped for ease of notation.

Up until now the noise has been kept entirely general. However in animaging system different types of noise are present. In its simplestform, the noise can be broken up into camera specific noise such as readout noise, dark current, etc. and photon noise or shot noise. Thevariance in Eq. (11) can thus be written as

$\begin{matrix}{{\sigma_{est}^{2} = {\frac{1}{N^{2}}{\sum\limits_{j = 0}^{N - 1}{\sum\limits_{k = 0}^{N - 1}{{\psi_{j}}^{- 2}\left\lbrack {{\sigma_{cam}^{2} + \sigma_{k}^{2}},{SN}} \right\rbrack}}}}},} & {{Equation}\mspace{14mu} (12)}\end{matrix}$

where σ² _(cam) is the variance of the camera noise (and does not dependon the index k since it is constant) and σ² _(k, SN) is the variance ofthe shot noise. This does depend on the index k since the shot noise canvary from image to image as it depends explicitly on the number ofphoton collected. In other words, imaging system 100 may change filter120 for different images, thereby inducing the number of photonscollected changes and affecting the shot noise. However for broadbandsources, this variation is small.

Eq. (12) may be simplified by assuming a broadband scene such that thenumber of photons collected from image to image is roughly constant.This allows us to take σ_(k,SN) ²≈ σ _(SN) ²w where the overbar denotessome average value. The noise estimate may be bounded by taking|ψ_(j)|=|ψ_(min)|. These simplifications give

$\begin{matrix}{{\sigma_{est}^{2} \lesssim {\frac{K_{2}}{n^{2}}\left\lbrack {\sigma_{cam}^{2} + {\overset{\_}{\sigma}}_{SN}^{2}} \right\rbrack}},} & {{Equation}\mspace{14mu} (13)}\end{matrix}$

where

$K = \frac{\psi_{\max}}{\psi_{\min}}$

is the spectral condition number using the fact that ψ_(max)=Σ_(i=0)^(N-1) c_(i)=n for a binary circulant matrix, where n is the number of1s in a row of the circulant matrix.

The variance term on the left had side of Eq. (13) is the variance ofthe signal estimate for a particular spectral band, while the terms onthe right hand side are variances of the image noise which contain noisefrom all wavelength bands, and thus depend on the kind of filter used.To compare the performance of this approach to the performance of asingle passband filter we will to convert the variance terms on theright hand side of Eq. (13) into terms that contain no explicitdependence on the type of filter used. The camera noise is alwaysconstant, and thus the noise is independent of the filter. However, theshot noise varies depending on how many photons reach the focal planearray, and that depends on how many spectral lines are open in filter120. Assuming simple relationship where the image shot noise is relatedto the total shot noise in the form

$\begin{matrix}{{{\overset{\_}{\sigma}}_{SN}^{2} \approx {\frac{n}{N}\sigma_{TSN}^{2}}},} & {{Equation}\mspace{14mu} (14)}\end{matrix}$

where σ² _(TSN) is the total shot noise if no filter was present and thefilter function S(λ)=1. This is the limiting case of a broadbandpanchromatic camera. Thus, Eq. (13) becomes

$\begin{matrix}{\sigma_{est}^{2} \lesssim {{\frac{K_{2}}{n_{2}}\left\lbrack {\sigma_{cam}^{2} + {\frac{n}{N}\sigma_{TSN}^{2}}} \right\rbrack}.}} & {{Equation}\mspace{14mu} (15)}\end{matrix}$

Compare this result to the noise one obtained using a simple passbandfilter. In that case the noise is

$\begin{matrix}{\sigma_{PB}^{2} = {\sigma_{cam}^{2} + {\frac{1}{N}{\sigma_{TSN}^{2}.}}}} & {{Equation}\mspace{14mu} (16)}\end{matrix}$

The gain G is defined as a ratio of the two SNR termsG=SNR_(est)/SNR_(PB)=σ_(PB)/σ_(est), where the second equality comesfrom the fact that the signal at a given band is equal in both thepassband and coded filter cases. The ratio G is then

$\begin{matrix}{{G = \frac{{n\left( {{N\; \chi^{2}} + 1} \right)}^{1/2}}{{K\left( {{N\; \chi^{2}} + n} \right)}^{1/2}}},} & {{Equation}\mspace{14mu} (17)}\end{matrix}$

where χ²=σ² _(cam)/σ² _(TSN) is the ratio of the camera to total shotnoise variance. The goal is to create filters such that G>1 whichimplies that imaging system 100 with tunable filter 120 and circulantQRS coded aperture 130 gives a higher SNR than a single passband filter.Here, G is inversely proportional to the condition number K as expected,matrices with a low condition number are desirable. However, it is alsodesirable for the ratio n/N to be as close to 1 as possible, to obtain ahighest possible SNR for each measurement. Unfortunately these twolimits cannot coincide since the condition number of a matrix typicallyscales as sparsity decreases (i.e., as the number of 0's in the matrixgoes down, the condition number goes up). To obtain an analytic result,that has the optimal SNR when χ>1, imaging system 100 implements aquadratic residue filter. The condition number K is known and given byK=√{square root over (2n)}, and n is related to N as n=(N+1)/2. Thus thegain G is given by

$\begin{matrix}{G = {\left\lbrack \frac{\left( {N + 1} \right)\left( {{N\; \chi^{2}} + 1} \right)}{2\left( {{2\; N\; \chi^{2}} + N + 1} \right)} \right\rbrack^{1/2}\overset{N1}{\rightarrow}{\sqrt{N}{\left( \frac{\chi^{2}}{{4\; \chi^{2}} + 1} \right)^{1/2}.}}}} & {{Equation}\mspace{14mu} (18)}\end{matrix}$

The final limit for large N shows that, for imaging system 100 whichimplements a quadratic residue filter with a large number of spectralbands, the standard deviation of the noise is reduced relative to thesingle passband case by a factor of √{square root over (N)}. Inaddition, this advantage increases as the signal becomes weaker, becauseG becomes larger for large values of χ.

FIG. 9 illustrates a graph of the gain ratio given in Eq. (18) for avariety of values of χ. FIG. 9 shows that the gain in dB is typicallygreater than 0 implying that imaging system 100 is an improvement over apassband approach. As seen in the graphic, the improvement is best forlarge values χ. Consider the gain improvement for large values of χ. Thestandard deviation is reduced in the multiple passband case by a factorof √{square root over (N)} relative to the single passband case becausemultiple images of the same spectral line are taken and averagedtogether in the process of taking the matrix inverse. Since the noise isadditive Gaussian noise, the standard √{square root over (N)} reductionin standard deviation is achieved. Shot noise is slightly more subtle,but the net payoff is the same.

Since more filter bands are open in imaging system 100, each image hasmore shot noise than if a single passband was used, but the signal doesnot increase, since the added signal and shot noise is at wavelength sdifferent from the signal of interest. This adds noise by a factor of√{square root over (N)} eliminating the advantage one gains by adding Nimages of the same spectral line together. However, as the number ofbands increases the shot noise per spectral band must decrease. Since afixed amount of light is divided up amongst more spectral bands for asingle passband filter, imaging system 100 still gains an SNRimprovement by a factor of √{square root over (N)} since the spectrum isdivided among fewer spectral bands.

Consider also a comparison of imaging system 100 which implements aquadratic residue spectral filter to that of a simple panchromaticcamera with no spectral information. Panchromatic imagers have the bestSNR characteristics of any camera, and are useful for remote sensingwhen low photon flux limits the ability to obtain spectral information.The shot noise variance in the panchromatic camera has already beendefined as σ² _(TSN), and the camera noise variance, σ² _(cam) will staythe same (assuming an identical camera is used). The approximation inEq. (14) relates the shot noise in a given spectral band to the totalshot noise available. Because the signal and shot noise are related byPoisson statistics, the same estimate t may be used to relate the totalsignal in a panchromatic imager to the signal at a given spectral linesuch that s_(pc)=(N/n) si. Here s_(pc) is the signal of a panchromaticimager, and s_(i) is the average signal at a spectral line λ_(i). Underthese conditions the SNR gain is

$\begin{matrix}{G = {{\sqrt{\frac{N + 1}{2\; N}}\left( \frac{\chi^{2} + 1}{{2\; N\; \chi^{2}} + N + 1} \right)^{1/2}}\overset{N1}{\rightarrow}{\frac{1}{\sqrt{2\; N}}{\left( \frac{\chi^{2} + 1}{{2\; \chi^{2}} + 1} \right)^{1/2}.}}}} & {{Equation}\mspace{14mu} (19)}\end{matrix}$

Thus we see that the panchromatic imager does indeed have better SNRthan imaging system 100 in this case. While imaging system 100 has alower SNR than a panchromatic imager, it provides a √{square root over(N)} improvement in the signal to noise over a single passband filter.Imaging system 100 thus provides the ability to do hyperspectral imagingin a relatively low light environment. Imaging system 100 also enablesuse of a high speed shutter to image transient events or to have veryhigh spectral resolution (many spectral bands).

Three Band Multiplexed Spectral Imager Embodiment

Consider, for the purposes of the discussion to follow, that imagingsystem 100 implements a three band multiplexed spectral imager. Asdescribed in greater detail below, processor 150 may further improve thegain by exploiting correlation properties of the noise and a prioriinformation about the scene being imaged. Further assume for thisexample that imaging system 100 or filter 120 implements or comprises aquadratic residue filter. To simplify the notation consider just onepixel of detector 140, and denote the measurement at this pixel bym_(i). To relate this to the imaging equation given in Eq. (5),m_(i)=I_(i) (u, v) and s_(i)=U″_(o)(u, v, λ_(i)) for a specific value of(u, v).

The measurement comprises the signals at the various spectral bandsdenoted by s_(i) and noise denoted by n_(i). Imaging system 100 takesthree measurements

m ₁ =s ₁ +s ₂ +n ₁ m ₂ =s ₂ +s ₃ +n ₂ m ₃ =s ₁ +s ₃ =n ₃  Equation (20):

The estimators for the signals are given in Eq. (10), and in the threeband case they simplify to

$\begin{matrix}{{e_{1} = {{\frac{1}{2}\left( {m_{1} - m_{2} + m_{3}} \right)} = {s_{1} + {\frac{1}{2}\left( {n_{1} - n_{2} + n_{3}} \right)}}}}{e_{2} = {{\frac{1}{2}\left( {m_{1} + m_{2} - m_{3}} \right)} = {s_{2} + {\frac{1}{2}\left( {n_{1} + n_{2} - n_{3}} \right)}}}}{e_{3} = {{\frac{1}{2}\left( {{- m_{1}} + m_{2} + m_{3}} \right)} = {s_{3} + {\frac{1}{2}\left( {{- n_{1}} + n_{2} + n_{3}} \right)}}}}} & {{Equation}\mspace{14mu} (21)}\end{matrix}$

The variance of each estimate in Eq. (21) is the same, and is given by

$\begin{matrix}{\sigma_{e}^{2} = {\frac{1}{4}{\left( {\sigma_{1}^{2} + \sigma_{2}^{2} + \sigma_{3}^{2}} \right).}}} & {{Equation}\mspace{14mu} (22)}\end{matrix}$

In the case where the variance of each measurement is the same (e.g.,when camera noise dominates shot noise) this simplifies to

$\sigma_{e}^{2} = {\frac{3}{4}{\sigma^{2}.}}$

This is the previous best case.

The noise terms of Eq. (21) are anti-correlated. The process of addingm₁+m₂+m₃ results in various noise terms cancelling each other because ofthis anti-correlation. Thus, if imaging system 100 performs measurementsin anti-correlated noise then multiple measurements of the signal willgive an SNR increase that is greater than the standard √{square rootover (N)} increase because of this noise-canceling effect. Previousimaging systems are unable to exploit this anti-correlation since thesignal is not known.

Imaging system 100 is also enhanced to use a priori information toexploit this anti-correlation. Even though the signal is unknown,imaging system 100 may use a priori information of the ratios betweenvarious signals. Assuming the ratios between signals are obtained, thesignals may be represented as s₂=α₂s₁ and s₃=α₃s₁. To find the optimalestimate such that the variance is minimized, define

$\begin{matrix}{{\hat{e_{1}} = {{{w_{1}e_{1}} + {w_{2}e_{2}} + {w_{3}e_{3}}} = {{\left( {w_{1} + {w_{2}\alpha_{2}} + {w_{3}\alpha_{3}}} \right)s_{1}} + {\frac{w_{1}}{2}\left( {n_{1} - n_{2} + n_{3}} \right)} + {\frac{w_{2}}{2}\left( {n_{1} + n_{2} - n_{3}} \right)} + {\frac{w_{3}}{2}\left( {{- n_{1}} + n_{2} + n_{3}} \right)}}}},} & {{Equation}\mspace{14mu} (23)}\end{matrix}$

where w_(i) is a weighting factor with the constraint thatw₁+w₂α₂+w₃α₃=1, so that

is in fact an estimator for s_(i). A similar definition and operationcould be performed for

or

but no generality is lost since we assume that we know the ratiosbetween the various signals that we seek to estimate.

The variance of

is given by

$\begin{matrix}{\sigma_{\hat{e}}^{2} = {{\frac{\sigma_{1}^{2}}{4}\left( {w_{1} + w_{2} - w_{3}} \right)} + {\frac{\sigma_{2}^{2}}{4}\left( {{- w_{1}} + w_{2} + w_{3}} \right)} + {\frac{\sigma_{3}^{2}}{4}{\left( {w_{1} - w_{2} + w_{3}} \right).}}}} & {{Equation}\mspace{14mu} (24)}\end{matrix}$

The variance of this estimator may be minimized subject to theconstraint that w₁+w₂α₂+w₃α₃=1. This results in a function

f(w ₁ ,w ₂ ,w ₃)=σ² _(ê)+λ(w ₁ +w ₂α₂ +w ₃α₃−1),  Equation (25):

where λ is the Lagrange multiplier. The minimization of this functionresults in a system of four equations and four unknowns (3 for the w'sand 1 for the constraint). The solutions are

$\begin{matrix}{{w_{1} = {\frac{1}{D}\left( {{\sigma_{1}^{2}\sigma_{2}^{2}} + {\sigma_{2}^{2}\sigma_{3}^{2}} + {\alpha_{2}\sigma_{2}^{2}\sigma_{3}^{2}} + {\alpha_{3}\sigma_{1}^{2}\sigma_{2}^{2}}} \right)}}{w_{2} = {\frac{1}{D}\left( {{\sigma_{2}^{2}\sigma_{3}^{2}} + {\alpha_{2}\left( {{\sigma_{1}^{2}\sigma_{3}^{2}} + {\sigma_{2}^{2}\sigma_{3}^{2}}} \right)} + {\alpha_{3}\sigma_{1}^{2}\sigma_{3}^{2}}} \right)}}{w_{1} = {\frac{1}{D}\left( {{\sigma_{1}^{2}\sigma_{2}^{2}} + {\alpha_{2}\sigma_{1}^{2}\sigma_{3}^{2}} + {\alpha_{3}\left( {{\sigma_{1}^{2}\sigma_{1}^{2}} + {\sigma_{1}^{2}\sigma_{3}^{2}}} \right)}} \right)}}} & {{Equation}\mspace{14mu} (26)}\end{matrix}$

where D is given by

D=(1+α₂)²σ² ₂σ² ₃+(α₂+α₃)²σ² ₁σ² ₃+(1+α₃)²σ² ₁σ² ₂,  Equation (27):

These solutions yield a variance for the estimator in Eq. (23) of

$\begin{matrix}{\sigma_{\hat{e}}^{2} = {\frac{\sigma_{1}^{2}\sigma_{2}^{2}\sigma_{3}^{2}}{D}.}} & {{Equation}\mspace{14mu} (28)}\end{matrix}$

If all variances are equal and if α₂=α₃=1, then Eq. (28) simplifies to

$\sigma_{\hat{e}}^{2} = {\frac{\sigma^{2}}{12}.}$

This is a dramatic improvement over the previous best case of

${\sigma_{\hat{e}}^{2} = {\frac{3}{4}\sigma^{2}}},$

and a significant improvement over what one would achieve by measuringthe signal 3 times.

FIG. 10 is a graph that illustrates the SNR gain when exploiting theanti-correlated nature of the noise. The horizontal axis is the ratior=α₂=α₃ and the vertical axis is the gain in dB calculated as 10 log G.The results show that using the estimator in Eq. (23) results in adramatic improvement in SNR over the estimator of Eq. (21).

Unfortunately in practice it is difficult to know the ratios α₂ and α₃exactly. In fact, if the estimates for the signals given in Eq. (21) areused to calculate the ratios, such that α₂=e₂/e₁ and α₃=e₃/e₁ then theseestimates combined with Eqs. (23) and (26) simply give

=e₁. In other words, the new estimator for the signal is no better thanour initial estimator. This makes sense since no additional informationhas been obtained in this instance, and the initial assumption that theratios have been obtained has been ignored.

To gain the dramatic improvement implied by Eq. (28), processor 150 mayobtain a better measure for the ratios. In practice, a scene has spatialcorrelations. Thus, the spectrum of a given pixel in a scene iscorrelated to the spectrum of its neighboring pixels. Imaging system 100may use this correlation information to generate a better estimate ofthe spectral ratios.

Assume that the spectral lines are related in the following way: s₂=α₂s₁+∈₂ and s₃=α₃ s₁+∈₃ where Σ_(i) is an error parameter that is includedif the estimate of α_(i) is not exactly correct. Further assume that∈_(i) is a Gaussian random variable with zero mean and variance σ² _(ε).Also, take the variance of ∈_(i) equal for all i for notationalsimplicity. Just as before, a general estimate of signal 1 may definedas

$\begin{matrix}{{\hat{e_{1}} = {{{w_{1}e_{1}} + {w_{2}e_{2}} + {w_{3}e_{3}}} = {{\left( {w_{1} + {w_{2}\alpha_{2}} + {w_{3}\alpha_{3}}} \right)s_{1}} + {\frac{w_{1}}{2}\left( {n_{1} - n_{2} + n_{3}} \right)} + {\frac{w_{2}}{2}\left( {{n_{1} + n_{2} - n_{3} + 2} \in_{21}} \right)} + {\frac{w_{3}}{2}\left( {{{- n_{1}} + n_{2} + {n_{3 +}2}} \in_{31}} \right)}}}},} & {{Equation}\mspace{14mu} (29)}\end{matrix}$

whose variance is

$\begin{matrix}{\sigma_{\hat{e}}^{2} = {{\frac{\sigma_{1}^{2}}{4}\left( {w_{1} + w_{2} - w_{3}} \right)^{2}} + {\frac{\sigma_{2}^{2}}{4}\left( {{- w_{1}} + w_{2} + w_{3}} \right)^{2}} + {\frac{\sigma_{3}^{2}}{4}\left( {w_{1} - w_{2} + w_{3}} \right)^{2}} + {w_{2}^{2}\sigma_{\varepsilon}^{2}} + {w_{3}^{2}{\sigma_{\varepsilon}^{2}.}}}} & {{Equation}\mspace{14mu} (30)}\end{matrix}$

To minimize Eq. (30) with respect to w₁, w₂, and w₃ yet simplify theanalysis, take σ² ₁=σ² ₂=σ² ₃=σ², and define the ratio χ=σ_(ε)/σ. Thevariance may now be minimized when

$\begin{matrix}{{w_{1} = {\frac{1}{D}\left\lbrack {{2\left( {2 + \alpha_{2} + \alpha_{3}} \right)} + {\left( {6 + \alpha_{2} + \alpha_{3}} \right)\chi^{2}} + {2\chi^{4}}} \right\rbrack}}{w_{2} = {\frac{1}{D}\left\lbrack {{2\left( {1 + {2\alpha_{2}} + \alpha_{3}} \right)} + {\left( {1 + {3\alpha_{2}}} \right)\chi^{2}}} \right\rbrack}}{w_{3} = {\frac{1}{D}\left\lbrack {{2\left( {1 + \alpha_{12} + {2\alpha_{3}}} \right)} + {\left( {1 + {3\alpha_{3}}} \right)\chi^{2}}} \right\rbrack}}} & {{Equation}\mspace{14mu} (31)}\end{matrix}$

where D is given by

D=2[(1+α₂)²+(α₂+α₃)²+(1+α₃)²]+[6+2(α₂+α₃)+3(α₂₂+α₂₃)]χ²+2χ⁴.  Equation(32):

The minimized variance in this case is now

$\begin{matrix}{\sigma_{\hat{e}}^{2} = {\frac{\sigma^{2}}{D^{2}}{\left\{ {{4\left\lbrack {\left( {1 + \alpha_{2}} \right)^{2} + \left( {\alpha_{2} + \alpha_{3}} \right)^{2} + \left( {1 + \alpha_{3}} \right)^{2}} \right\rbrack} + {\left\lbrack {24 + {8\left( {\alpha_{2} + \alpha_{3}} \right)} + {12\left( {\alpha_{2}^{2} + \alpha_{3}^{2}} \right)}} \right\rbrack \chi^{2}} + {2\left( {15 + \alpha_{2} + {3\alpha_{2}^{2}} + \alpha_{3} + {3\alpha_{3}^{2}} - {3\alpha_{2}\alpha_{3}}} \right)\chi^{4}} + {16\chi^{6}} + {3\chi^{8}}} \right\}.}}} & {{Equation}\mspace{14mu} (33)}\end{matrix}$

These solutions may be examined in two limits. First when χ>>1, thisimplies that σ_(∈)>>σ. In other words, the error on the estimate of theratios α₂ and α₃ is very high. In this limit we see that w₁=1, w₂=0,w_(S)=0, w₃=0, and σ_(ê)=3/4σ².

FIG. 11 illustrates a graphic of the SNR gain as a function of χ. Thehorizontal axis is the ratio χ=σ_(ε)/σ and the vertical axis is the gainin dB calculated as 10 log G. The results show that using the estimatorin Eq. (23) results in a dramatic improvement in SNR over the estimatorof Eq. (21). The gain fails off as χ increases. For very large values ofχ that gain is still significant. Even for χ=20 there is a gain ofroughly 5 dB.

The three band example discussed above may be generalized for usefulapplications of imaging system 100. The same notation is retained forsimplicity. The coded N band measurement is given in Eq. (5). It may berewritten in single pixel notation as

m _(i)=Σ_(j=0) ^(N-1) q _(ij) s _(j) +n _(i).  Equation (34):

The estimator for the signal is given in Eq. (10), and is rewritten insimplified notation as

e ₁ =s ₁+Σ_(k=0) ^(N-1) r _(k-1) n _(k),  Equation (35):

where

$\begin{matrix}{r_{k - l} = {\frac{1}{N}{\sum\limits_{j = 0}^{N - 1}{\psi_{j}^{- 1}{^{{- 2}\pi \; {{({k - 1})}}{j/N}}.}}}}} & {{Equation}\mspace{14mu} (36)}\end{matrix}$

For the case of a quadratic residue filter, Eq. (36) reduces to

$r_{k - l} = {\frac{2}{N + 1}\left( {{2q_{k - 1}} - 1} \right)}$

as shown in Eq. (8).

Just as with the three band example, we now assume that the ratiosbetween the spectral bands are known. Thus we write all bands in termsof band 0 as

s _(i)=α_(i) s ₀+∈_(i)  Equation (37):

where π_(i)=s_(i)/s_(o) is the ratio between signal s_(i) and s_(o) and∈_(i) accounts for any error in the estimate in the ratio. New estimators_(o) is created by combining all of the estimators of the various bandswith a weighting factor such that

=Σ_(l=0) ^(N-1) w ₁ e ₁=Σ_(l=0) ^(N-1)(w ₁α_(i) s _(o)+∈_(i)+Σ_(k=0)^(N-1) r _(k-1) n _(k)),  Equation (38):

subject to the constraint that Σ_(l=0) ^(N-1) W₁α₁=1. The variance ofthis new estimator is given by

π² _(ê)=Σ_(l=0) ^(N-1) =w ² ₁σ_(∈,L) ²+Σ_(k=0) ^(N-1)(Σ_(l=0) ^(N-1) w_(l) r _(k-l))²σ² _(k),  Equation (39):

where σ_(∈,l) ₁ ²=

∈² ₁

is the expectation value of the square error in the ratio estimate.

To minimize Eq. (39) with the constraint given above, the function maybe defined as

f(w _(l),λ)≡σ_(ê) ²+λ(Σ_(l=0) ^(N-1) w _(l)α_(l)−1),  Equation (40):

which may be minimized with respect to w_(i) and λ. Taking the gradientwith respect to the w_(l)'s and λ yields the following N+1 with N+1unknowns

2w _(l)σ_(∈,l) ²+2Σ_(k=0) ^(N-1)Σ_(j=0) ^(N-1) w _(j) r _(k-j) r_(k-l)σ_(k) ²+λα_(l)=0

Σ_(t=0) ^(N-1) w _(l)α_(l)=1.  Equation (41):

Eq. (41) may be solved for the weights w1 so that an estimator iscreated given Eq. (38) that has the noise minimized. To solve Eq. (41)in general requires an N+1×N+1 matrix inverse operation for each pixelon the focal plane array. This problem quickly becomes too timeconsuming for more than a few spectral bands and pixels. However a fewsimplifying assumptions may be taken which allow the matrix inverseoperation to be done very quickly. The sacrifice is that the true noiseminimum will not be found, but numerical tests have shown that the errorincrease is minimal while the speed increase is dramatic.

The first simplifying assumption regards the noise terms σ_(k) ²=σ_(cam)²+σ² _(k, SN). Assume that the camera noise is much greater than theshot noise such that σ² _(k)=σ_(cam) ²≡σ² is independent of theparameter k. Assume further that error term σ_(∈,l) ² may be ignored.This term accounts for the uncertainty in the exact knowledge of thespectral ratios. Ignoring this term shows the worst case scenario (i.e.,simply the original estimator), which occurs when no a priori knowledgeof the ratios is used. Thus the error will generally not get worse thanthe original estimator. Numerical simulations have shown that ignoringthis term does cause a slight decrease in the optimization routine, butagain this loss is minimal compared to the speed gained when ignoringthis term. With these assumptions Eq. (41) may be rewritten as

Σ_(k=0) ^(N-1)Σ_(j=0) ^(N-1) w _(j) r _(k-j) r _(k-l)+λ′α_(l)=0 Σ_(l=0)^(N-1) =w _(l)α_(l)=1.  Equation(42):

A new Lagrange multiplier is defined as λ′=λ/2σ² for simplicity. Thereason for these simplifications becomes apparent when the equation iswritten in matrix form

$\begin{matrix}{{\begin{pmatrix}{r_{k}r_{k}} & {r_{k - 1}r_{k}} & \ldots & {r_{k - {({N - 1})}}r_{k}} & \alpha_{0} \\{r_{k}r_{k - 1}} & {r_{k - 1}r_{k - 1}} & \ldots & {r_{k - {({N - 1})}}r_{k - 1}} & \alpha_{1} \\\vdots & \vdots & \ddots & \vdots & \vdots \\{r_{k}r_{k - {({N - 1})}}} & {r_{k - 1}r_{k - {({N - 1})}}} & \ldots & {r_{k - {({N - 1})}}r_{k - {({N - 1})}}} & \alpha_{N - 1} \\\alpha_{0} & \alpha_{1} & \ldots & \alpha_{N - 1} & \alpha_{0}\end{pmatrix}\begin{pmatrix}w_{0} \\w_{1} \\\vdots \\w_{N - 1} \\\lambda^{\prime}\end{pmatrix}} = {\begin{pmatrix}0 \\0 \\\vdots \\0 \\1\end{pmatrix}.}} & {{Equation}\mspace{14mu} (43)}\end{matrix}$

Eq. (43) uses summation notation such that any two variablesside-by-side that contain the same index are summed over that index.Thus r_(k-j)r_(k)=Σ_(k=0) ^(N-1) r_(k-j)r_(k). Negative indices wraparound as do indices great that N−1.

The four blocks of Eq. (43) are now defined. The N×N upper left block is

$\begin{matrix}{R = \begin{pmatrix}{r_{k}r_{k}} & {r_{k - 1}r_{k}} & \ldots & {r_{k - {({N - 1})}}r_{k}} \\{r_{k}r_{k - 1}} & {r_{k - 1}r_{k - 1}} & \ldots & {r_{k - {({N - 1})}}r_{k - 1}} \\\vdots & \vdots & \ddots & \vdots \\{r_{k}r_{k - {({N - 1})}}} & {r_{k - 1}r_{k - {({N - 1})}}} & \ldots & {r_{k - {({N - 1})}}r_{k - {({N - 1})}}}\end{pmatrix}} & {{Equation}\mspace{14mu} (44)}\end{matrix}$

The vector a is

$\begin{matrix}{{a = \begin{pmatrix}\alpha_{0} \\\alpha_{1} \\\vdots \\\alpha_{N - 1}\end{pmatrix}},} & {{Equation}\mspace{14mu} (45)}\end{matrix}$

along with a^(T) its transpose, and the vector w is

$\begin{matrix}{w = {\begin{pmatrix}w_{0} \\w_{1} \\\vdots \\w_{N - 1} \\\lambda^{\prime}\end{pmatrix}.}} & {{Equation}\mspace{14mu} (46)}\end{matrix}$

These definitions allow Eq. (43) to be rewritten as

$\begin{matrix}{{\begin{pmatrix}R & a \\a^{T} & 0\end{pmatrix}\begin{pmatrix}w \\\lambda^{\prime}\end{pmatrix}} = {\begin{pmatrix}0 \\1\end{pmatrix}.}} & {{Equation}\mspace{14mu} (47)}\end{matrix}$

Processor 150 may remove the 0 in the bottom right term to invert Eq.(47). Processor 150 may do this by adding the row directly above thelast row since it is equal to 0. This yields an equation in the form

$\begin{matrix}{{{\begin{pmatrix}R & a \\b^{T} & \alpha_{N - 1}\end{pmatrix}\begin{pmatrix}w \\{\lambda^{\prime} + w_{N - 1}}\end{pmatrix}} = \begin{pmatrix}0 \\1\end{pmatrix}},} & {{Equation}\mspace{14mu} (48)}\end{matrix}$

where b^(T) is the vector a^(T) with the row directly above it added toit. Processor 150 may invert the matrix in Block form, which gives thesolution for the w vector as

$\begin{matrix}{w = {{- \frac{1}{c}}R^{- 1}a}} & {{Equation}\mspace{14mu} (49)}\end{matrix}$

where c=α_(N-1)−b^(T) R^(×1) a.

The solution for the weighting terms w_(i) now requires an N×N inverseof the matrix R along with a matrix multiplication step. The reason forthe speed increase is that the matrix R is constant over the entirefocal plane, thus the inverse only needs to be calculated once. Only thevectors a and b^(T) change for each pixel. In addition the matrix R iscirculant, and thus its inverse can be calculated using a Fouriertransform. Thus, processor 150 may invert the N×N block using a FastFourier Transform (FFT). We have now reduced our problem fromcalculating an N+1×N+1 matrix inverse for every single pixel on ourfocal plane to calculate a length N Fourier transform once, and doing asimple matrix multiplication for each pixel on the focal plane. Imagingsystem 100 is therefore enhanced to provide a dramatic speed increase tothe noise minimization algorithm.

As described above, imaging system 100 may use neighboring pixels tobetter estimate these ratios, and obtain large SNR enhancement. Thisrelies on the assumption that the spectrum varies slowly over a fewpixels on the focal plane. For natural scenes this should typically bethe case. If this assumption is not valid other methods may also existto estimate the ratios such as predefined spectral catalogs or physicalmodels.

FIG. 12 illustrates an exemplary sequence of images 1201-1203 taken withthe same filter which is circularly shifted for each image. Consider asimulated scene which contains four circular sources each with uniformintensity across the source, but with unique spectral signatures. FIG.13 illustrates an exemplary panchromatic sequence of images taken of ascene which contains four circular sources 1310, 1320, 1330, and 1340each with a unique spectral profile 1301-1303. The noise-freepanchromatic images correspond to an output one may expect using theexemplary filters of FIG. 13.

FIG. 14 illustrates the actual spectrum for each source of FIG. 13 in anexemplary embodiment. The upper left plot 1410 of FIG. 14 correspondswith the upper left source 1302 of FIG. 13 and so on. FIG. 15illustrates a high resolution reconstruction of the spectrum from thepanchromatic images of FIG. 13 in an exemplary embodiment. Here,processor 150 used a length 31 quadratic residue and expanded it to alength 310 array in the simulation thus providing a 310 spectral bandresolution.

FIG. 16 illustrates a reconstruction of the spectrum from thepanchromatic images of FIG. 13 in an exemplary embodiment. Here, imageprocesser 150 has reconstructed the spectrum from only 31 images. ThoughFIG. 16 shows resolution degradation since there are only 31 spectralbands the peaks are nonetheless clearly visible. It will be appreciatedthat this process is merely exemplary for what happens when using afewer number of images to do the reconstruction and thus does notrepresent a best possible result.

FIG. 17 illustrates a noisy image in an exemplary embodiment. Thesources 1710-1740 of the scene 1701 are captured in an image whichincludes shot noise and camera noise. Assume for this example, that theshot noise is taken to be a Gaussian random variable with varianceproportional to the intensity of the image, and the camera noise is alsoa Gaussian random variable but with constant variance from image toimage.

Consider reconstruction of the spectral datacube of a simulated scene,in the presence of the noise as in FIG. 17 using prior techniques (e.g.,single passband filters) as well as using the processes described hereinwith respect to imaging system 100. FIG. 18 illustrates a comparison ofspectrums for a pixel in the noisy image of FIG. 17 in an exemplaryembodiment. Plot 1801 shows the actual spectrum of a pixel of the upperleft source 1710 in the image of FIG. 17. Plot 1802 shows areconstructed spectrum of the pixel using 107 single passband filters.Plot 1803 shows a reconstructed spectrum using 107 quadratic residuesequences in an exemplary embodiment (e.g., reconstructed with theprocesses of imaging system 100). For this example, assume a filterresolution with a length 107 quadratic residue sequence with a 1 nmGaussian passband. In this simulation, the peak signal of the cameranoise ratio (ignoring shot noise) is 20 if no filter is used and allphotons of all wavelength s are captured.

Plot 1803 shows a dramatic improvement when using the quadratic residuefilter. Using Eq. (18), and the noise parameters used in thissimulation, the estimated SNR obtained when using a quadratic residuecoded filter should be at least 3.7 times that of the single passbandcase or roughly a 5.5 dB improvement in SNR. Compare this estimate tothe standard deviation of the actual signal minus the passband signal,and the actual signal minus the signal using imaging system 100. The 5.5dB gain estimate agrees quantitatively with the ratio of these errormeasures, and is qualitatively consistent with the improvement seen inFIG. 18

For some scenes, such as natural scenes, it may be assumed that thespectrum varies slowly from pixel to pixel. Processor 150 may thus beconfigured to estimate the ratios using a nearest neighbor average. Inthis example, processor 150 uses five nearest neighbors of a pixel. FIG.19 illustrates a comparison of spectrums for a pixel in the noisy imageof FIG. 17 in another exemplary embodiment. Here, plot 1901 shows theactual spectrum of a pixel of the upper left source 1710 in the image ofFIG. 17. And plot 1902 shows a reconstructed spectrum using 107quadratic residue sequences in an exemplary embodiment (e.g.,reconstructed with the processes of imaging system 100) and with thenoise minimization algorithm applied.

Experimental Results

FIG. 20 illustrates images obtained from simulated measurements of ahyperspectral datacube in an exemplary embodiment. The simulations ofthis example used actual hyperspectral image data obtained from the NASAMASTER program, a joint development involving the Airborne SensorFacility at the Ames Research Center, the Jet Propulsion Laboratory andthe EROS Data Center that provides imaging data online for publicavailability. Plot 2010 shows an RGB image of the hyperspectral “truth”datacube. Plots 2020 shows the simulated effect of measuring thedatacube one band at a time using a passband filter, and plot 2030 showsthe simulated effect of measuring multiple bands using a length 47quadratic residue filter. Plot 2040 shows the effect of performing thenoise minimization algorithm on the simulated noisy datacube obtainedfrom using the quadratic reside filter. Thus, plot 2040 shows a blackand white version of the false color RGB image obtained after performingthe noise minimization. FIG. 20 shows an improvement in the imagequality over that shown in plots 1320 and 1330, although the gain may bedifficult to see in print form.

Here, the five nearest neighbors of a pixel were used to estimate thespectral ratios. The relative error for each datacube compared to thetruth was calculated in each case. The SNR gain from changing from thepassband to QR filter was roughly a factor of 2 or about 3 dB. The SNRgain one obtains when performing the noise minimization is a factor or 7better than the QR measurement or about 8.5 dB improvement.

FIG. 21 illustrates a comparison of the spectrum for the central pixelof FIG. 20. The plots of FIG. 21 highlight the dramatic SNR gain. Plot2110 shows the truth point spectrum, plot 2120 shows the simulatedpassband point spectrum measurement, plot 2130 shows the simulated QRpoint spectrum measurement, and plot 2140 shows the simulated QRmeasurement with noise reduction algorithm. Each spectrum becomesprogressively better, resulting in the final measurement which has verylittle noise corruption.

The invention can also take the form of an entirely hardware embodiment,an entirely software embodiment or an embodiment containing bothhardware and software elements. In one embodiment, the invention isimplemented in software, which includes but is not limited to firmware,resident software, microcode, etc. FIG. 22 illustrates a computingsystem 2200 in which a computer readable medium 2206 may provideinstructions for performing any of the methods disclosed herein.

Furthermore, the invention can take the form of a computer programproduct accessible from the computer readable medium 2206 providingprogram code for use by or in connection with a computer or anyinstruction execution system. For the purposes of this description, thecomputer readable medium 2206 can be any apparatus that can tangiblystore the program for use by or in connection with the instructionexecution system, apparatus, or device, including the computer system2200.

The medium 2206 can be any tangible electronic, magnetic, optical,electromagnetic, infrared, or semiconductor system (or apparatus ordevice). Examples of a computer readable medium 2206 include asemiconductor or solid state memory, magnetic tape, a removable computerdiskette, a random access memory (RAM), a read-only memory (ROM), arigid magnetic disk and an optical disk. Some examples of optical disksinclude compact disk-read only memory (CD-ROM), compact disk-read/write(CD-R/W) and DVD.

The computing system 2200, suitable for storing and/or executing programcode, can include one or more processors 2202 coupled directly orindirectly to memory 2208 through a system bus 2210. The memory 2208 caninclude local memory employed during actual execution of the programcode, bulk storage, and cache memories which provide temporary storageof at least some program code in order to reduce the number of timescode is retrieved from bulk storage during execution. Input/output orI/O devices 2204 (including but not limited to keyboards, displays,pointing devices, etc.) can be coupled to the system either directly orthrough intervening I/O controllers. Network adapters may also becoupled to the system to enable the computing system 2200 to becomecoupled to other data processing systems, such as through host systemsinterfaces 2212, or remote printers or storage devices throughintervening private or public networks. Modems, cable modem and Ethernetcards are just a few of the currently available types of networkadapters.

What is claimed is:
 1. A spectral imaging system, comprising: an opticalelement configured to receive electromagnetic energy of atwo-dimensional scene; a filter configured to provide a plurality ofspectral filter profiles, and to transmit multiple spectral wavebands ofthe electromagnetic energy substantially simultaneously through at leastone of the plurality of spectral profiles; a detector configured tomeasure intensities of the multiple spectral wavebands transmittedthrough the at least one spectral profile; and a processor configured togenerate a spectral image of the scene based on the measuredintensities.
 2. The system of claim 1, wherein: the spectral filterprofiles each comprise a static multiple-waveband filter; the detectorcomprises a single sensor element of a two-dimensional array of sensorelements; and each static multiple-waveband filter transmits multiplespectral wavebands onto individual sensor elements of thetwo-dimensional array.
 3. The system of claim 1, further comprising:multiple detectors arranged at separate image planes from one another;wherein the spectral filter profiles comprise static multiple-wavebandfilters that transmit different portions of the electromagnetic energyfrom one another; and wherein each static multiple-waveband filtertransmits a portion of the electromagnetic energy on a different one ofthe multiple detectors.
 4. The system of claim 3, wherein: at least oneof the static multiple-waveband filters transmits the electromagneticenergy in a different direction than another one of the staticmultiple-waveband filters.
 5. The system of claim 1, further comprising:a coded aperture configured to multiplex the transmitted spectralwavebands; wherein the detector comprises a single sensor element of atwo-dimensional array of sensor elements; and wherein the single sensorelement is configured to measure two or more of the multiplexed spectralwavebands substantially simultaneously.
 6. The system of claim 5,wherein: the coded aperture comprises a matrix in the form of atwo-dimensional pattern of sections; and each section blocks or reflectssome or all of the impinging electromagnetic energy to transform thetransmitted spectral wavebands on the detector.
 7. The system of claim1, wherein: the filter is configured to sequentially modify one or moreof the spectral profiles in response to an input, and to spectrallyencode the multiple spectral wavebands in a matrix pattern.
 8. Thesystem of claim 7, wherein: the matrix pattern includes invertibleproperties; and the processor is configured to decode the multiplexedspectral wavebands on the single sensor element into individual spectralcomponents by inverting the matrix based on the invertible properties,and to generate the spectral image based with the individual spectralcomponents.
 9. The system of claim 8, wherein: the matrix is a circulantmatrix that includes a quadratic residue sequence.
 10. The system ofclaim 1, further comprising: the processor is configured to reduceerrors in the spectral image based on a ratio between two or more of themeasured intensities;
 11. The system of claim 11, wherein: the processoris configured to obtain the ratio based on spatial correlations in thespectral image.
 12. The system of claim 11, wherein: the processor isconfigured to obtain the ratio based on a-prior information of thescene.
 13. The system of claim 1, wherein: the electromagnetic energyincludes a red waveband, a blue waveband, and a green waveband; and thefilter comprises: a first spectral profile configured to transmit thered waveband and the blue waveband, and to attenuate the green waveband;a second spectral profile configured to transmit the red waveband andthe green waveband, and to attenuate the blue waveband; and a thirdspectral profile configured to transmit the green waveband and the bluewaveband, and to attenuate the red waveband.
 14. A method for imaging ascene, the method comprising: receiving electromagnetic energy of atwo-dimensional scene; transmitting, with a filter, multiple spectralwavebands of the electromagnetic energy substantially simultaneouslythrough at least one of a plurality of spectral profiles of the filter;measuring, with a detector, intensities of the multiple spectralwavebands; and generating a spectral image of the scene based on themeasured intensities.
 15. The method of claim 14, further comprising:measuring, with a single sensor element in a two-dimensional array ofthe detector, two or more of the transmitted wavebands substantiallysimultaneously.
 16. The method of claim 14, wherein: the spectral filterprofiles comprise static multiple-waveband filters that transmitdifferent portions of the electromagnetic energy from one another; andeach static multiple-waveband filter transmits a portion of theelectromagnetic energy on a different detector.
 17. The method of claim13, further comprising: sequentially modifying, with the filter, one ormore of the spectral profiles in response to an input; spectrallyencoding the multiple spectral wavebands in a matrix pattern.
 18. Themethod of claim 14, wherein: the matrix is a circulant matrix thatincludes a quadratic residue sequence.
 19. The method of claim 18,wherein: the matrix pattern includes invertible properties; and themethod further comprises: decoding the multiplexed spectral wavebands onthe single sensor element into individual spectral components byinverting the matrix based on the invertible properties; and generatingthe spectral image based with the individual spectral components. 20.The method of claim 16, further comprising: reducing errors in thespectral image based on a ratio between two or more of the measuredintensities.